/*
  File:    NecrosisLevelSets.cpp
  Purpose: Routines for maintaining necrotic and quiescent level set functions. 

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cstdio>
#include <cstdlib>
#include <iostream>

using namespace std;

#include <cstdlib>
#include <cmath>

using namespace std;

#include "../Matrix.h"
#include "../ImprovedMath.h"

#include "Contour.h"
#include "DiscreteReinit.h"
#include "reinit.h"

int DiscreteSign( double input )
{
 extern double _EPS;

 if( input > _EPS )
 { return 1; }
 if( input < -_EPS )
 { return -1; }
 return 0;
}

void UpdateNecrosisLevelset( void )
{
 extern double _EPS;
 
 extern int Xnodes;
 extern int Ynodes;
 extern double dx;
 extern double dy;
 extern Matrix X;
 extern Matrix Y;
 
 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 extern Matrix Oxygen;
 
 extern double NecrosisOxygenLevel;
 
 if( PhiNecrotic.TellRows() < 2 )
 { PhiNecrotic = ones(Xnodes,Ynodes); }
 
 if( Oxygen.TellRows() < 2 )
 { Oxygen = ones(Xnodes,Ynodes); }
 
 // find horizontal points on the interface. 

 Matrix Temporary = zeros(Xnodes,Ynodes);
 
 int i,j;
 for( j=0 ; j < Oxygen.TellCols() ; j++ )
 {
  for( i=0 ; i < Oxygen.TellRows() ; i++ )
  {
   bool done = false;
   if( *Phi(i,j) <= _EPS && *PhiNecrotic(i,j) <= _EPS )
   { *Temporary(i,j) = -1.0; done = true; }
   if( *Phi(i,j) <= _EPS && *Oxygen(i,j) - NecrosisOxygenLevel <= _EPS )
   { *Temporary(i,j) = -1.0; done = true; }
   if( !done )
   { *Temporary(i,j) = 1.0; done = true; }
  }
 }
 
 Contour ZeroContour;
 DataPoint* pCurrentDataPoint = ZeroContour.pFirstDataPoint;
 
 for( j=0 ; j < Oxygen.TellCols() ; j++ )
 {
  for( i=0 ; i < Oxygen.TellRows()-1 ; i++ )
  {
   if( *Temporary(i,j) > 0 && *Temporary(i+1,j) < 0 )
   {
    double phi_center = *Phi(i,j);
	double phi_right  = *Phi(i+1,j);
	
	double phi_necrotic_center = *PhiNecrotic(i,j);
	double phi_necrotic_right  = *PhiNecrotic(i+1,j);
   
    double sig_center = *Oxygen(i,j) - NecrosisOxygenLevel;
	double sig_right  = *Oxygen(i+1,j) - NecrosisOxygenLevel;
   
    double ThetaPhi = fabs(phi_center) / ( fabs(phi_center) + fabs(phi_right) );
	if( phi_center <= _EPS && phi_right <= _EPS )
	{ ThetaPhi = 0; }
	
	double ThetaPhiN = fabs(phi_necrotic_center) / 
	( fabs(phi_necrotic_center) + fabs(phi_necrotic_right) );
	if( phi_necrotic_center > _EPS && phi_necrotic_right > _EPS )
	{ ThetaPhiN = 2; }
	if( phi_necrotic_center <= _EPS && phi_necrotic_right <= _EPS )
	{ ThetaPhiN = 0; }
	if( ThetaPhiN < ThetaPhi )
	{ ThetaPhiN = ThetaPhi; }
	
	double ThetaPhiSigma = fabs(sig_center) / ( fabs(sig_center) + fabs(sig_right) );
	if( sig_center > _EPS && sig_right > _EPS )
	{ ThetaPhiSigma = 2; }
	if( sig_center <= _EPS && sig_right <= _EPS )
	{ ThetaPhiSigma = 0; }
	if( ThetaPhiSigma < ThetaPhi )
	{ ThetaPhiSigma = ThetaPhi; }
	
	double Theta = ThetaPhiN; 
	if( ThetaPhiSigma < Theta )
	{ Theta = ThetaPhiSigma; }
	
	DataPoint* pNewDataPoint;
	pNewDataPoint = new DataPoint;
	pNewDataPoint->x = *X(i) + Theta*dx;
	pNewDataPoint->y = *Y(j);
	
	if( pCurrentDataPoint )
	{ pCurrentDataPoint->pNextDataPoint = pNewDataPoint; }
	else
	{ ZeroContour.pFirstDataPoint = pNewDataPoint; }
	pCurrentDataPoint = pNewDataPoint;
   }

   if( *Temporary(i,j) < 0 && *Temporary(i+1,j) > 0 )
   {
    double phi_center = *Phi(i,j);
	double phi_right  = *Phi(i+1,j);
	
	double phi_necrotic_center = *PhiNecrotic(i,j);
	double phi_necrotic_right  = *PhiNecrotic(i+1,j);
   
    double sig_center = *Oxygen(i,j) - NecrosisOxygenLevel;
	double sig_right  = *Oxygen(i+1,j) - NecrosisOxygenLevel;
   
    double ThetaPhi = fabs(phi_center) / ( fabs(phi_center) + fabs(phi_right) );
	if( phi_center <= _EPS && phi_right <= _EPS )
	{ ThetaPhi = 1; }
	
	double ThetaPhiN = fabs( phi_necrotic_center ) / 
	( fabs(phi_necrotic_center) + fabs(phi_necrotic_right) );
	if( phi_necrotic_center > _EPS && phi_necrotic_right > _EPS )
	{ ThetaPhiN = -2; }
	if( phi_necrotic_center <= _EPS && phi_necrotic_right <= _EPS )
	{ ThetaPhiN = 1; }
	if( ThetaPhi < ThetaPhiN )
	{ ThetaPhiN = ThetaPhi; }
	
	double ThetaPhiSigma = fabs( sig_center ) / ( fabs(sig_center) + fabs(sig_right) );
	if( sig_center > _EPS && sig_right > _EPS )
	{ ThetaPhiSigma = -2; }
	if( sig_center <= _EPS && sig_right <= _EPS )
	{ ThetaPhiSigma = 1; }
	if( ThetaPhi < ThetaPhiSigma )
	{ ThetaPhiSigma = ThetaPhi; }
	
	double Theta = ThetaPhiN; 
	if( ThetaPhiSigma > Theta )
	{ Theta = ThetaPhiSigma; }
	
	DataPoint* pNewDataPoint;
	pNewDataPoint = new DataPoint;
	pNewDataPoint->x = *X(i) + Theta*dx;
	pNewDataPoint->y = *Y(j);
	
	if( pCurrentDataPoint )
	{ pCurrentDataPoint->pNextDataPoint = pNewDataPoint; }
	else
	{ ZeroContour.pFirstDataPoint = pNewDataPoint; }
	pCurrentDataPoint = pNewDataPoint;
   }
   
  }
 }
 
 for( i=0 ; i < Oxygen.TellRows() ; i++ )
 {
  for( j=0 ; j < Oxygen.TellCols()-1 ; j++ )
  {
   if( *Temporary(i,j) > 0 && *Temporary(i,j+1) < 0 )
   {
    double phi_center = *Phi(i,j);
	double phi_right  = *Phi(i,j+1);
	
	double phi_necrotic_center = *PhiNecrotic(i,j);
	double phi_necrotic_right  = *PhiNecrotic(i,j+1);
   
    double sig_center = *Oxygen(i,j) - NecrosisOxygenLevel;
	double sig_right  = *Oxygen(i,j+1) - NecrosisOxygenLevel;
   
    double ThetaPhi = fabs(phi_center) / ( fabs(phi_center) + fabs(phi_right) );
	if( phi_center <= _EPS && phi_right <= _EPS )
	{ ThetaPhi = 0; }
	
	double ThetaPhiN = fabs( phi_necrotic_center ) / 
	( fabs(phi_necrotic_center) + fabs(phi_necrotic_right) );
	if( phi_necrotic_center > _EPS && phi_necrotic_right > _EPS )
	{ ThetaPhiN = 2; }
	if( phi_necrotic_center <= _EPS && phi_necrotic_right <= _EPS )
	{ ThetaPhiN = 0; }
	if( ThetaPhiN < ThetaPhi )
	{ ThetaPhiN = ThetaPhi; }
	
	double ThetaPhiSigma = fabs( sig_center ) / ( fabs(sig_center) + fabs(sig_right) );
	if( sig_center > _EPS && sig_right > _EPS )
	{ ThetaPhiSigma = 2; }
	if( sig_center <= _EPS && sig_right <= _EPS )
	{ ThetaPhiSigma = 0; }
	if( ThetaPhiSigma < ThetaPhi )
	{ ThetaPhiSigma = ThetaPhi; }
	
	double Theta = ThetaPhiN; 
	if( ThetaPhiSigma < Theta )
	{ Theta = ThetaPhiSigma; }
	
	DataPoint* pNewDataPoint;
	pNewDataPoint = new DataPoint;
	pNewDataPoint->x = *X(i);
	pNewDataPoint->y = *Y(j) + Theta*dy;
	
	if( pCurrentDataPoint )
	{ pCurrentDataPoint->pNextDataPoint = pNewDataPoint; }
	else
	{ ZeroContour.pFirstDataPoint = pNewDataPoint; }
	pCurrentDataPoint = pNewDataPoint;
   }

   if( *Temporary(i,j) < 0 && *Temporary(i,j+1) > 0 )
   {
    double phi_center = *Phi(i,j);
	double phi_right  = *Phi(i,j+1);
	
	double phi_necrotic_center = *PhiNecrotic(i,j);
	double phi_necrotic_right  = *PhiNecrotic(i,j+1);
   
    double sig_center = *Oxygen(i,j) - NecrosisOxygenLevel;
	double sig_right  = *Oxygen(i,j+1) - NecrosisOxygenLevel;
   
    double ThetaPhi = fabs(phi_center) / ( fabs(phi_center) + fabs(phi_right) );
	if( phi_center <= _EPS && phi_right <= _EPS )
	{ ThetaPhi = 1; }
	
	double ThetaPhiN = fabs( phi_necrotic_center ) / 
	( fabs(phi_necrotic_center) + fabs(phi_necrotic_right) );
	if( phi_necrotic_center > _EPS && phi_necrotic_right > _EPS )
	{ ThetaPhiN = -2; }
	if( phi_necrotic_center <= _EPS && phi_necrotic_right <= _EPS )
	{ ThetaPhiN = 1; }
	if( ThetaPhi < ThetaPhiN )
	{ ThetaPhiN = ThetaPhi; }
	
	double ThetaPhiSigma = fabs( sig_center ) / ( fabs(sig_center) + fabs(sig_right) );
	if( sig_center > _EPS && sig_right > _EPS )
	{ ThetaPhiSigma = -2; }
	if( sig_center <= _EPS && sig_right <= _EPS )
	{ ThetaPhiSigma = 1; }
	if( ThetaPhi < ThetaPhiSigma )
	{ ThetaPhiSigma = ThetaPhi; }
	
	double Theta = ThetaPhiN; 
	if( ThetaPhiSigma > Theta )
	{ Theta = ThetaPhiSigma; }
	
	DataPoint* pNewDataPoint;
	pNewDataPoint = new DataPoint;
	pNewDataPoint->x = *X(i);
	pNewDataPoint->y = *Y(j) + Theta*dy;
	
	if( pCurrentDataPoint )
	{ pCurrentDataPoint->pNextDataPoint = pNewDataPoint; }
	else
	{ ZeroContour.pFirstDataPoint = pNewDataPoint; }
	pCurrentDataPoint = pNewDataPoint;
   }
   
  }
 }
  
// PhiNecrotic = Reinit( Temporary ); //  DiscreteReinit( Temporary , ZeroContour );
 PhiNecrotic = DiscreteReinit( Temporary , ZeroContour );
// PhiNecrotic = Reinit( PhiNecrotic );
// PhiNecrotic = Reinit( PhiNecrotic );
 
 return;
}

void UpdateQuiescentLevelset( void )
{
 extern double _EPS;
 
 extern Matrix X;
 extern Matrix Y;
 
 extern int Xnodes;
 extern int Ynodes;
 extern double dx;
 extern double dy;
 
 extern Matrix Phi;
 extern Matrix PhiQuiescent;
 extern Matrix Oxygen;
 
 extern double QuiescenceOxygenLevel;
 
 if( PhiQuiescent.TellRows() < 2 )
 { PhiQuiescent = ones(Xnodes,Ynodes); }
 
 if( Oxygen.TellRows() < 2 )
 { Oxygen = ones(Xnodes,Ynodes); }
 
 // find horizontal points on the interface. 

 Matrix Temporary = zeros(Xnodes,Ynodes);
 
 int i,j;
 for( j=0 ; j < Oxygen.TellCols() ; j++ )
 {
  for( i=0 ; i < Oxygen.TellRows() ; i++ )
  {
   bool done = false;
   if( *Phi(i,j) <= _EPS && *PhiQuiescent(i,j) <= _EPS )
   { *Temporary(i,j) = -1.0; done = true; }
   if( *Phi(i,j) <= _EPS && *Oxygen(i,j) - QuiescenceOxygenLevel <= _EPS )
   { *Temporary(i,j) = -1.0; done = true; }
   if( !done )
   { *Temporary(i,j) = 1.0; done = true; }
  }
 }
 
 Contour ZeroContour;
 DataPoint* pCurrentDataPoint = ZeroContour.pFirstDataPoint;
 
 for( j=0 ; j < Oxygen.TellCols() ; j++ )
 {
  for( i=0 ; i < Oxygen.TellRows()-1 ; i++ )
  {
   if( *Temporary(i,j) > 0 && *Temporary(i+1,j) < 0 )
   {
    double phi_center = *Phi(i,j);
	double phi_right  = *Phi(i+1,j);
	
	double phi_quiescent_center = *PhiQuiescent(i,j);
	double phi_quiescent_right  = *PhiQuiescent(i+1,j);
   
    double sig_center = *Oxygen(i,j) - QuiescenceOxygenLevel;
	double sig_right  = *Oxygen(i+1,j) - QuiescenceOxygenLevel;
   
    double ThetaPhi = fabs(phi_center) / ( fabs(phi_center) + fabs(phi_right) );
	if( phi_center <= _EPS && phi_right <= _EPS )
	{ ThetaPhi = 0; }
	
	double ThetaPhiQ = fabs(phi_quiescent_center) / 
	( fabs(phi_quiescent_center) + fabs(phi_quiescent_right) );
	if( phi_quiescent_center > _EPS && phi_quiescent_right > _EPS )
	{ ThetaPhiQ = 2; }
	if( phi_quiescent_center <= _EPS && phi_quiescent_right <= _EPS )
	{ ThetaPhiQ = 0; }
	if( ThetaPhiQ < ThetaPhi )
	{ ThetaPhiQ = ThetaPhi; }
	
	double ThetaPhiSigma = fabs(sig_center) / ( fabs(sig_center) + fabs(sig_right) );
	if( sig_center > _EPS && sig_right > _EPS )
	{ ThetaPhiSigma = 2; }
	if( sig_center <= _EPS && sig_right <= _EPS )
	{ ThetaPhiSigma = 0; }
	if( ThetaPhiSigma < ThetaPhi )
	{ ThetaPhiSigma = ThetaPhi; }
	
	double Theta = ThetaPhiQ; 
	if( ThetaPhiSigma < Theta )
	{ Theta = ThetaPhiSigma; }
	
	DataPoint* pNewDataPoint;
	pNewDataPoint = new DataPoint;
	pNewDataPoint->x = *X(i) + Theta*dx;
	pNewDataPoint->y = *Y(j);
	
	if( pCurrentDataPoint )
	{ pCurrentDataPoint->pNextDataPoint = pNewDataPoint; }
	else
	{ ZeroContour.pFirstDataPoint = pNewDataPoint; }
	pCurrentDataPoint = pNewDataPoint;
   }

   if( *Temporary(i,j) < 0 && *Temporary(i+1,j) > 0 )
   {
    double phi_center = *Phi(i,j);
	double phi_right  = *Phi(i+1,j);
	
	double phi_quiescent_center = *PhiQuiescent(i,j);
	double phi_quiescent_right  = *PhiQuiescent(i+1,j);
   
    double sig_center = *Oxygen(i,j) - QuiescenceOxygenLevel;
	double sig_right  = *Oxygen(i+1,j) - QuiescenceOxygenLevel;
   
    double ThetaPhi = fabs(phi_center) / ( fabs(phi_center) + fabs(phi_right) );
	if( phi_center <= _EPS && phi_right <= _EPS )
	{ ThetaPhi = 1; }
	
	double ThetaPhiQ = fabs( phi_quiescent_center ) / 
	( fabs(phi_quiescent_center) + fabs(phi_quiescent_right) );
	if( phi_quiescent_center > _EPS && phi_quiescent_right > _EPS )
	{ ThetaPhiQ = -2; }
	if( phi_quiescent_center <= _EPS && phi_quiescent_right <= _EPS )
	{ ThetaPhiQ = 1; }
	if( ThetaPhi < ThetaPhiQ )
	{ ThetaPhiQ = ThetaPhi; }
	
	double ThetaPhiSigma = fabs( sig_center ) / ( fabs(sig_center) + fabs(sig_right) );
	if( sig_center > _EPS && sig_right > _EPS )
	{ ThetaPhiSigma = -2; }
	if( sig_center <= _EPS && sig_right <= _EPS )
	{ ThetaPhiSigma = 1; }
	if( ThetaPhi < ThetaPhiSigma )
	{ ThetaPhiSigma = ThetaPhi; }
	
	double Theta = ThetaPhiQ; 
	if( ThetaPhiSigma > Theta )
	{ Theta = ThetaPhiSigma; }
	
	DataPoint* pNewDataPoint;
	pNewDataPoint = new DataPoint;
	pNewDataPoint->x = *X(i) + Theta*dx;
	pNewDataPoint->y = *Y(j);
	
	if( pCurrentDataPoint )
	{ pCurrentDataPoint->pNextDataPoint = pNewDataPoint; }
	else
	{ ZeroContour.pFirstDataPoint = pNewDataPoint; }
	pCurrentDataPoint = pNewDataPoint;
   }
   
  }
 }
 
 for( i=0 ; i < Oxygen.TellRows() ; i++ )
 {
  for( j=0 ; j < Oxygen.TellCols()-1 ; j++ )
  {
   if( *Temporary(i,j) > 0 && *Temporary(i,j+1) < 0 )
   {
    double phi_center = *Phi(i,j);
	double phi_right  = *Phi(i,j+1);
	
	double phi_quiescent_center = *PhiQuiescent(i,j);
	double phi_quiescent_right  = *PhiQuiescent(i,j+1);
   
    double sig_center = *Oxygen(i,j) - QuiescenceOxygenLevel;
	double sig_right  = *Oxygen(i,j+1) - QuiescenceOxygenLevel;
   
    double ThetaPhi = fabs(phi_center) / ( fabs(phi_center) + fabs(phi_right) );
	if( phi_center <= _EPS && phi_right <= _EPS )
	{ ThetaPhi = 0; }
	
	double ThetaPhiQ = fabs( phi_quiescent_center ) / 
	( fabs(phi_quiescent_center) + fabs(phi_quiescent_right) );
	if( phi_quiescent_center > _EPS && phi_quiescent_right > _EPS )
	{ ThetaPhiQ = 2; }
	if( phi_quiescent_center <= _EPS && phi_quiescent_right <= _EPS )
	{ ThetaPhiQ = 0; }
	if( ThetaPhiQ < ThetaPhi )
	{ ThetaPhiQ = ThetaPhi; }
	
	double ThetaPhiSigma = fabs( sig_center ) / ( fabs(sig_center) + fabs(sig_right) );
	if( sig_center > _EPS && sig_right > _EPS )
	{ ThetaPhiSigma = 2; }
	if( sig_center <= _EPS && sig_right <= _EPS )
	{ ThetaPhiSigma = 0; }
	if( ThetaPhiSigma < ThetaPhi )
	{ ThetaPhiSigma = ThetaPhi; }
	
	double Theta = ThetaPhiQ; 
	if( ThetaPhiSigma < Theta )
	{ Theta = ThetaPhiSigma; }
	
	DataPoint* pNewDataPoint;
	pNewDataPoint = new DataPoint;
	pNewDataPoint->x = *X(i);
	pNewDataPoint->y = *Y(j) + Theta*dy;
	
	if( pCurrentDataPoint )
	{ pCurrentDataPoint->pNextDataPoint = pNewDataPoint; }
	else
	{ ZeroContour.pFirstDataPoint = pNewDataPoint; }
	pCurrentDataPoint = pNewDataPoint;
   }

   if( *Temporary(i,j) < 0 && *Temporary(i,j+1) > 0 )
   {
    double phi_center = *Phi(i,j);
	double phi_right  = *Phi(i,j+1);
	
	double phi_quiescent_center = *PhiQuiescent(i,j);
	double phi_quiescent_right  = *PhiQuiescent(i,j+1);
   
    double sig_center = *Oxygen(i,j) - QuiescenceOxygenLevel;
	double sig_right  = *Oxygen(i,j+1) - QuiescenceOxygenLevel;
   
    double ThetaPhi = fabs(phi_center) / ( fabs(phi_center) + fabs(phi_right) );
	if( phi_center <= _EPS && phi_right <= _EPS )
	{ ThetaPhi = 1; }
	
	double ThetaPhiQ = fabs( phi_quiescent_center ) / 
	( fabs(phi_quiescent_center) + fabs(phi_quiescent_right) );
	if( phi_quiescent_center > _EPS && phi_quiescent_right > _EPS )
	{ ThetaPhiQ = -2; }
	if( phi_quiescent_center <= _EPS && phi_quiescent_right <= _EPS )
	{ ThetaPhiQ = 1; }
	if( ThetaPhi < ThetaPhiQ )
	{ ThetaPhiQ = ThetaPhi; }
	
	double ThetaPhiSigma = fabs( sig_center ) / ( fabs(sig_center) + fabs(sig_right) );
	if( sig_center > _EPS && sig_right > _EPS )
	{ ThetaPhiSigma = -2; }
	if( sig_center <= _EPS && sig_right <= _EPS )
	{ ThetaPhiSigma = 1; }
	if( ThetaPhi < ThetaPhiSigma )
	{ ThetaPhiSigma = ThetaPhi; }
	
	double Theta = ThetaPhiQ; 
	if( ThetaPhiSigma > Theta )
	{ Theta = ThetaPhiSigma; }
	
	DataPoint* pNewDataPoint;
	pNewDataPoint = new DataPoint;
	pNewDataPoint->x = *X(i);
	pNewDataPoint->y = *Y(j) + Theta*dy;
	
	if( pCurrentDataPoint )
	{ pCurrentDataPoint->pNextDataPoint = pNewDataPoint; }
	else
	{ ZeroContour.pFirstDataPoint = pNewDataPoint; }
	pCurrentDataPoint = pNewDataPoint;
   }
   
  }
 }
 
// PhiQuiescent = Reinit( Temporary ); // DiscreteReinit( Temporary , ZeroContour );
 PhiQuiescent = DiscreteReinit( Temporary , ZeroContour );
// PhiQuiescent = Reinit( PhiQuiescent );
// PhiQuiescent = Reinit( PhiQuiescent );
  
 return;
}

Matrix NecroticLevelsetAdvectionOperator( void )
{
 extern double _EPS;
 
 extern int Xnodes;
 extern int Ynodes;
 extern double dx;
 extern double dy;
 
 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 
 extern Matrix Pressure;
 

 Matrix Output = zeros(Xnodes,Ynodes);
 
 int i,j;
 for( j=0 ; j<Ynodes ; j++ )
 {
  for( i=0 ; i < Xnodes ; i++ )
  {
   if( fabs( *PhiNecrotic(i,j) ) <= 5*dx )
   {
    double Px = 0;
	double Py = 0;
	
	bool done = false;
	if( DiscreteSign(*Phi(i-1,j))*DiscreteSign(*Phi(i+1,j)) < 1 )
	{ Px = (*Pressure(i+1,j)-*Pressure(i-1,j))/(2.0*dx); done = true; }
	if( !done && *Phi(i-1,j) <= _EPS )
	{ Px = (*Pressure(i,j) - *Pressure(i-1,j))/dx ; done = true; }
	if( !done && *Phi(i+1,j) <= _EPS )
	{ Px = (*Pressure(i+1,j) - *Pressure(i,j))/dx ; done = true; }
	
	done = false;
	if( *Phi(i,j-1) <= _EPS && *Phi(i,j+1) <= _EPS )
	{ Py = (*Pressure(i,j+1)-*Pressure(i,j-1))/(2.0*dy); done = true; }
	if( !done && *Phi(i,j-1) <= _EPS )
	{ Py = (*Pressure(i,j) - *Pressure(i,j-1))/dy ; done = true; }
	if( !done && *Phi(i,j+1) <= _EPS )
	{ Py = (*Pressure(i,j+1) - *Pressure(i,j))/dy ; done = true; }
   
   }
  
  
   if( fabs( *PhiNecrotic(i,j) ) <= 5*dx )
   {
    *Output(i,j) =  (  (*Pressure(i+1,j)-*Pressure(i-1,j))*
	                   (*PhiNecrotic(i+1,j)-*PhiNecrotic(i-1,j))
					  /(4*square(dx)) 
					  +(*Pressure(i,j+1)-*Pressure(i,j-1))*
					   (*PhiNecrotic(i,j+1)-*PhiNecrotic(i,j-1))
					  /(4*square(dy)) );
   }   
  }
 }

 return Output;
}

